#####################################
## "Backyard politics in Foreign Aid"
## William Christiansen
## Tobias Heinrich
## Timothy Peterson
#####################################

## File that generates tables of
## all estimated models


## Load outputs
###############
load("output/est_out_YM.Rdata")
load("output/est_out_YT.Rdata")
load("output/est_out_MT.Rdata")
load("output/data_d1.Rdata")
load("output/data_d2.Rdata")


## Table for Y | T
##################
mod1 <- lm(Y1 ~ 1 + T1 + T2, data=subset(d1, Change == "Increase"))
mod2 <- lm(Y2 ~ 1 + T1 + T2, data=subset(d2, Change == "Increase"))
mod3 <- lm(Y1 ~ 1 + T1 + T2, data=subset(d1, Change == "Cut"))
mod4 <- lm(Y2 ~ 1 + T1 + T2, data=subset(d2, Change == "Cut"))

out_YT_tmp <- ddply(.data= melt(out_YT), .variables=c("DV", "Scenario", "Coeffcient"),
                    .fun=function(x) data.frame(m=mean(x$value),
                                                l=quantile(x$value, .025), u=quantile(x$value, .975)))
out_YT_tmp <- subset(out_YT_tmp, m != 0)
out_YT_tmp$Coeffcient2 <- pretty_var_name(i=as.character(out_YT_tmp$Coeffcient))

tab <- stargazer(mod1, mod2, mod3, mod4,
                 coef=list(subset(out_YT_tmp, Scenario == "Increase" & DV == "Rating")$m,
                           subset(out_YT_tmp, Scenario == "Increase" & DV == "Thermometer")$m,
                           subset(out_YT_tmp, Scenario == "Cut" & DV == "Rating")$m,
                           subset(out_YT_tmp, Scenario == "Cut" & DV == "Thermometer")$m),
                 ci=TRUE, digits=1, style="io", model.numbers=FALSE,
                 no.space = TRUE,
                 dep.var.labels=c("Project", "Senator", "Project", "Senator"),
                 column.labels=c("Increase foreign aid", "Cut foreign aid"), column.separate=c(2, 2), report="vcs",
                 font.size="footnotesize", table.placement="ht!", df=FALSE, omit.stat=c("ll", "adj.rsq", "f", "rsq"),
                 ci.custom=list(subset(out_YT_tmp, Scenario == "Increase" & DV == "Rating")[, c("l", "u")],
                                subset(out_YT_tmp, Scenario == "Increase" & DV == "Thermometer")[, c("l", "u")],
                                subset(out_YT_tmp, Scenario == "Cut" & DV == "Rating")[, c("l", "u")],
                                subset(out_YT_tmp, Scenario == "Cut" & DV == "Thermometer")[, c("l", "u")]),
                 covariate.labels=c("Local project, no Senator", "Local project, Senator", "Intercept"))
tab[10] <- paste0("\\\\[-1.8ex] & \\textbf{", paste0(rep(c("Project}", "Senator}"), 2), collapse=" & \\textbf{"), " \\\\ ")
tab <- c(tab[4], tab[7:22], 
         "\\hline \\\\[-1.8ex] ", "\\end{tabular} ",  
         "\\caption{\\textbf{Results for treatment effects.} Each single number presents the mean estimate; underneath is the 95\\% confidence interval. All models were bootstrapped; the Senator-models used the cluster-adjusted bootstrap.}\\label{Tab:YT}",
         "\\end{table}")
write(x=tab, file="output/R-Table-YT.tex")


## Tables for M | T, X
######################
for(i in 1:2)
{
  mod1 <- glm(formula=formula(paste0("M1 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  mod2 <- glm(formula=formula(paste0("M2 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  mod3 <- glm(formula=formula(paste0("M3 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  mod4 <- glm(formula=formula(paste0("M4 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  mod5 <- glm(formula=formula(paste0("M5 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  mod6 <- glm(formula=formula(paste0("M6 ~ 1", form_main, form_covs)), data=subset(d1, Change == ifelse(i == 1, "Increase", "Cut")), family=binomial())
  
  out_MT_tmp <- melt(out_MT[, , , "Expanded", ])
  out_MT_tmp <- ddply(.data=out_MT_tmp, .variables=c("M", "Scenario", "Coeffcient"),
                      .fun=function(x) data.frame(m=mean(x$value),
                                                  l=quantile(x$value, .025), u=quantile(x$value, .975)))
  out_MT_tmp$Coeffcient2 <- pretty_var_name(i=as.character(out_MT_tmp$Coeffcient))
  
  tab <- stargazer(mod1, mod2, mod3, mod4, mod5, mod6, 
                   coef=list(subset(out_MT_tmp, M == "M1" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m,
                             subset(out_MT_tmp, M == "M2" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m,
                             subset(out_MT_tmp, M == "M3" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m,
                             subset(out_MT_tmp, M == "M4" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m,
                             subset(out_MT_tmp, M == "M5" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m,
                             subset(out_MT_tmp, M == "M6" & Scenario == ifelse(i == 1, "Increase", "Cut"))$m),
                   ci=TRUE, digits=1, style="io", model.numbers=FALSE,
                   no.space = TRUE,
                   dep.var.labels=paste0("M", 1:6),
                   column.labels=NULL, report="vcs",
                   font.size="tiny", table.placement="ht!", df=FALSE, omit.stat=c("ll", "adj.rsq", "f", "rsq"),
                   ci.custom=list(subset(out_MT_tmp, M == "M1" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")],
                                  subset(out_MT_tmp, M == "M2" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")],
                                  subset(out_MT_tmp, M == "M3" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")],
                                  subset(out_MT_tmp, M == "M4" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")],
                                  subset(out_MT_tmp, M == "M5" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")],
                                  subset(out_MT_tmp, M == "M6" & Scenario == ifelse(i == 1, "Increase", "Cut"))[, c("l", "u")]),
                   covariate.labels=c("Local project, no Senator", "Local project, Senator", 
                                      "Gender", "Age", "Costs", "Education, high", "Education, low",  
                                      "Ideology", "Housing value", "Housing price increase", "Transnational",
                                      "'US greatest country'", "Intercept"))
  tab[10] <- paste0("\\\\[-1.8ex] & \\textbf{", paste0(c("Friends, family", "State", "US economy", "US security", 
                                                         "Elites", "Moral"), collapse="} & \\textbf{"), "} \\\\ ")
  tab <- c(tab[4], tab[7:51],
           "\\hline \\\\[-1.8ex] ", "\\end{tabular} ",
           paste0("\\caption{\\textbf{Results for mediators; aid ", ifelse(i == 1, "increase", "cut"), " cases only.} Coefficients and 95\\% confidence intervals from Bernoulli-probit models, 
                  regressing each mediator on a full covariate set.}\\label{Tab:MT-", ifelse(i == 1, "Increase", "Cut"), "}"),
           "\\end{table}")
  
  write(x=tab, file=paste0("output/R-Table-MT-", ifelse(i == 1, "Increase", "Cut"), ".tex"))
}  


## Table for Y | M, T, X
########################
this_form <- paste0("~ 1 +", paste0("M", 1:6, collapse="+"), form_main, form_covs)
mod1 <- lm(formula=formula(paste0("Y1", this_form)), data=subset(d1, Change == "Increase"))
mod2 <- lm(formula=formula(paste0("Y2", this_form)), data=subset(d2, Change == "Increase"))
mod3 <- lm(formula=formula(paste0("Y1", this_form)), data=subset(d1, Change == "Cut"))
mod4 <- lm(formula=formula(paste0("Y2", this_form)), data=subset(d2, Change == "Cut"))

out_YM_tmp <- ddply(.data=melt(out_YM), .variables=c("DV", "Scenario", "Coeffcient"),
                    .fun=function(x) data.frame(m=mean(x$value),
                                                l=quantile(x$value, .025), u=quantile(x$value, .975)))
out_YM_tmp$Coeffcient2 <- pretty_var_name(i=as.character(out_YM_tmp$Coeffcient))

tab <- stargazer(mod1, mod2, mod3, mod4,
                 coef=list(subset(out_YM_tmp, Scenario == "Increase" & DV == "Rating")$m,
                           subset(out_YM_tmp, Scenario == "Increase" & DV == "Thermometer")$m,
                           subset(out_YM_tmp, Scenario == "Cut" & DV == "Rating")$m,
                           subset(out_YM_tmp, Scenario == "Cut" & DV == "Thermometer")$m),
                 ci=TRUE, digits=1, style="io", model.numbers=FALSE,
                 no.space = TRUE,
                 dep.var.labels=c("Project", "Senator", "Project", "Senator"),
                 column.labels=c("Increase foreign aid", "Cut foreign aid"), column.separate=c(2, 2), report="vcs",
                 font.size="scriptsize", table.placement="ht!", df=FALSE, omit.stat=c("ll", "adj.rsq", "f", "rsq"),
                 ci.custom=list(subset(out_YM_tmp, Scenario == "Increase" & DV == "Rating")[, c("l", "u")],
                                subset(out_YM_tmp, Scenario == "Increase" & DV == "Thermometer")[, c("l", "u")],
                                subset(out_YM_tmp, Scenario == "Cut" & DV == "Rating")[, c("l", "u")],
                                subset(out_YM_tmp, Scenario == "Cut" & DV == "Thermometer")[, c("l", "u")]),
                 covariate.labels=c(paste0("M:", c("Friends, family", "State", "US economy", "US security",
                                                   "Elites", "Moral")), "Local project, no Senator", "Local project, Senator", 
                                    "Gender", "Age", "Costs", "Education, high", "Education, low", 
                                    "Ideology", "Housing value", "Housing price increase", "Transnational",
                                    "'US greatest country'", "Intercept"))
tab[10] <- paste0("\\\\[-1.8ex] & \\textbf{", paste0(rep(c("Project}", "Senator}"), 2), collapse=" & \\textbf{"), " \\\\ ")
tab <- c(tab[4], tab[7:70], 
          "\\hline \\\\[-1.8ex] ", "\\end{tabular} ", 
          "\\caption{\\textbf{Results for treatment effects using full covariate set (including mediators).} Each single number presents the mean estimate; underneath is the 95\\% confidence interval. All models were bootstrapped; the Senator-models used the cluster-adjusted bootstrap.}\\label{Tab:YMT}",
         "\\end{table} ")
write(x=tab, file="output/R-Table-YM.tex")
